Generalized Deam-Edwards Approach to the 
Statistical Mechanics of Randomly Crosslinked Systems 

Xiangjun Xing^, Bing-Sui Lu^, Fangfu Ye^, and Paul M. Goldbart^ 

^Department of Physics and Institute of Natural Sciences, 
Shanghai Jiao Tang University, Shanghai, China and 
'^School of Physics, Georgia Institute of Technology, 
837 State Street, Atlanta, GA 30332-0430, U.S.A. 
(Dated: May 6, 2013) 

We address the statistical mechanics of randomly and permanently crosslinked networks. We 
develop a theoretical framework (vulcanization theory) which can be used to systematically analyze 
the correlation between the statistical properties of random networks and their histories of forma- 
tion. Generalizing the original idea of Deam and Edwards, we consider an instantaneous crosslinking 
process, where all crosslinkers (modeled as Gaussian springs) are introduced randomly at once in 
an equilibrium liquid state, referred to as the preparation state. The probability that two func- 
tional sites are crosslinked by a spring exponentially decreases with their distance squared. After 
formally averaging over network connectivity, we obtained an effective theory with all degrees of 
freedom replicated 1 + n times. Two thermodynamic ensembles, the preparation ensemble and the 
measurement ensemble, naturally appear in this theory. The former describes the thermodynamic 
fluctuations in the state of preparation, while the latter describes the thermodynamic fluctuations in 
the state of measurement. We classify various correlation functions and discuss their physical signif- 
icances. In particular, the memory correlation functions characterize how the properties of networks 
depend on their history of formation, and are the hallmark properties of all randomly crosslinked 
materials. We clarify the essential difference between our approach and that of Deam-Edwards, 
discuss the saddle-point order parameters and its physical significance. Finally we also discuss the 
connection between saddle-point approximation of vulcanization theory, and the classical theory of 
rubber elasticity as well as the neo-classical theory of nematic elastomers. 

PACS numbers: 61.30.Vx,61.30.-v,61.43.-j 

I. INTRODUCTION 

Randomly crosslinked systems — including rubbers, gels, liquid crystalline elastomers, cytoskeleton networks, and 
many other related systems — arise widely in nature and are also extensively manufactured by scientists and tech- 
nologists . Their physical properties are dominated by heterogeneous, system-spanning random networks. A proper 
specification of this random network structure therefore constitutes the first step towards obtaining a more complete 
understanding of these complex materials. 

In a typical protocol for making randomly crosslinked materials, one starts from a liquid system under particular 
physical conditions (e.g., of temperature, pressure, solvents, etc.), and randomly introduces crosslinkers that link 
together nearby particles or polymers. The macroscopic liquid state prior to crosslinking is defined as the preparation 
state. A realistic crosslinking process necessarily takes some finite interval of time, and is therefore inevitably accom- 
panied by some irreversible responses of the system. Inclusion of these irreversible aspects would make the theoretical 
modeling of the random crosslinking process extremely complicated. To simplify the analysis and nevertheless capture 
the essential physics of random crosslinking, we shall therefore follow the original idea of Deam and Edwards [I] by 
considering an idealized, instantaneous crosslinking scheme in which all crosslinkers are introduced into the system at 
one particular instant in time. If a sufficient number of crosslinkers is introduced, a system-spanning random network 
emerges, and endows the material with a non-vanishing shear modulus. A continuous transition such as this, from 
a liquid to a random solid, is usually called gelation — or vulcanization in settings of pre-existing polymers. Gelation 
and vulcanization are quite different from usual thermodynamic phase transitions, because they are irreversible and 
the resulting systems acquire intrinsic spatial heterogeneity. The statistical theory for these transitions is therefore 
necessarily more complicated. 

It is well known that the physical properties of randomly crosslinked materials depend both on the state of mea- 
surement and on the state of preparation (i.e, the history of formation). This is, of course, a feature common to all 
heterogeneous materials, and has been widely appreciated by experimentalists. It was also explicitly pointed out by 
Deam and Edwards [T] and by de Gennes |^ as early as the 1970's. It therefore follows that the statistical physics of 
randomly crosslinked materials involves two statistical ensembles: the preparation ensemble and the measurement en- 
semble. This necessitates a substantial generalization of the standard Gibbs ensemble theory for equilibrium statistical 
mechanics. 
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To substantiate this seemingly extreme statement, let us discuss the remarkably distinct properties exhibited by 
nematic elastomers that have been crosslinked under distinct conditions. Nematic elastomers [3] are crosslinked 
nematic polymer melts, and their liquid crystalline ordering is strongly coupled to their network elasticity. Two 
crosslinking protocols have been studied recently by Urayama et al. In one, the system is first crosslinked in the high- 
temperature, isotropical phase of liquid crystallinity, and is then brought to a lower-temperature, nematic phase. The 
resulting system (called an isotropic-genesis polydomain (PD) nematic elastomer, or I-PNE) exhibits a polydomain 
nematic state, with a typical nematic domain size of a few microns. Upon stretching, the nematic domains gradually 
rotate and deform, and eventually, beyond certain strain value, they align along the external stress so as to form a 
monodomain (MD) state. The stress response during the PD-MD transition window is extremely soft, being orders 
of magnitude lower than that of typical elastomeric materials. 

In the other protocol, the system is first quenched into the nematic phase, so that the nematic order is well developed 
locally, but is still frustrated by nematic defects at larger scales (typically beyond 100 /im). Crosslinking carried out 
in this state freezes in the patterns of nematic defects. The resulting system (known as a nematic-genesis polydomain 
nematic elastomer, or N-PNE) also exhibits a stress-driven PD-MD transition. The observed stress is, however, within 
the range of typical elastomer materials, and much higher than that of I-PNEs. The striking differences between I- 
PNEs and N-PNEs are entirely due to their distinct preparation states. Evidently, a theoretical understanding of these 
differences demands a statistical theory that involves both the preparation ensemble and the measurement ensemble. 

A generic protocol for preparing and measuring a randomly crosslinked system has three steps, each statistical in 
nature. First, there is a preparation state. We shall denote averaging over the thermal fiuctuations in this state by 
square brackets [^]p , where the subscript p stands for "preparation'.' Second, there is an instantaneous crosslinking 
stage, at which crosslinkers are randomly introduced into the system. Averages over realizations of crosslinkers are 
denoted by an overbar: A . Finally, there is a measurement state. Averages over fluctuations in this state are 
denoted by angle brackets {A)^ , where the subscript m stands for "measurement'.' Note that these three types of 
average have a causal ordering in time: Any quantity averaged in a measurement state must also be averaged over 
crosslinker realizations and over the preparation state, as the latter two ensembles constitute the history of formation. 
By contrast, quantities averaged in the preparation ensemble are not to be averaged over crosslinking, because the 
network does not exist in the preparation state. 

From the above example of nematic elastomers, it can be inferred that the structure of a polymer network in 
a randomly crosslinked system depends not only on the randomness inherent in the crosslinking process but also 
on the thermal fluctuations in the state of preparation. Two macroscopically identical crosslinking processes using 
the identical protocol yield two distinct networks which are only statistically similar. More importantly, two such 
crosslinking processes carried out under different preparation states yield statistically different network structures, 
and lead to different properties of the systems, even under the same conditions in the measurement state. It therefore 
takes both the crosslinking methods and the preparation state to have a full statistical description of the structure of 
the random network in these materials. 

In view of the three steps just mentioned, as well as their statistical nature, there are four types of physical quan- 
tities that can be measured in a randomly crosslinked system: 

(1) Thermodynamic averages in the preparation state: [^]p, [^S]p etc., where the subscripts p stands for 
"preparation'.' 

(2) Thermodynamic averages in the measurement state: (^)m 7 {AB)^ etc. Note that these are to be averaged 
over both their crosslinking randomness and the thermal fluctuations in the preparation state. 



(3) Glassy correlations, which characterize sample-to-sample fiuctuations in the measurement state: {A)i-a{B) 
These quantities are important because the systems are random. 



(4) Memory correlations: A (-B),n , which characterize cross- correlations between the preparation and measurement 

states. Memory correlations therefore characterize the dependence of system properties on their histories of formation. 

The purpose of this Paper is to discuss various basic aspects of vulcanization theory, which provides a general 
methodology for the systematic calculation of all four of these types of diagnostic quantity. The remainder of the 
Paper is organized as follows. In Sec. [IT] we introduce a toy model of randomly linked anisotropic particles, as a 
starting point for vulcanization theory. In Sec. Ill we discuss a generalized Deam-Edwards distribution of network 
structure, and average the free energy of the network over this distribution. At this stage, we shall see a key theme, 
viz., that both the preparation and the measurement ensemble appear naturally, coupled to one another, in the 
effective replicated Hamiltonian. In Sec. |IV| we discuss in detail various physical quantities that can be calculated 
using the replicated partition function. We also discuss the significance of the zeroth replica and the connection 
between replica limit and causality. We also address the distinction between the original Deam-Edwards distribution 
and our generalization of it. In Sec. |V] we discuss the physical meaning of the vulcanization order parameter. In 
Sec.lVTlwe study the mean-field approximation, and from it we derive the classical theory of rubber elasticity, as well 



FIG. 1: Collection of particles, schematically depicted as roughly spherical blobs, containing functional sites labeled i and 
having spatial locations Ci. Some number of crosslinkers (here, Gaussian springs) are introduced into the system; each of them 
can link a pair of functional sites. The detailed interactions between the particles do not play a significant role in the present 
discussion. 

as the neoclassical theory elasticity, appropriate for nematic elastomers. Finally in Sec. |VII[ we draw the concluding 
remarks by briefly reviewing three basic results about the heterogeneous nature of randomly crosslinked materials, 
and by indicating some possible future directions. 

II. MODEL OF LINKED PARTICLES 

To make our discussion concrete, we consider a liquid of particles, each of which may or may not be anisotropic. 
The particles may, e.g., be entire polymers that we regard in a coarse-grained sense. Under appropriate conditions, 
the particles may participate in some kind of collective ordering, such as liquid crystalline alignment. Each particle 
has certain number of functional sites at which the crosslinkers can act. Before crosslinking, the spatial locations of 
all the functional sites of all the particles, which we index via i, are denoted by c^. We assume that the positions of 
these sites cP = {c°, i — 1, . . . , iV} completely specify the microscopic state of the liquid (i.e., there are no relevant 
degrees of freedom beyond these). The liquid Hamiltonian (normalized by temperature) is denoted by H^^^{c^). The 
corresponding partition function for the liquid system (in the preparation state) is then given (with the shorthand 
notation Dc^ = 11^=1 '^^1.) 

4, = /^c°e-^".(^"'. (1) 



TABLE I: Two equivalent specifications of the network structure shown in Fig. [T] 



X ^(1,9) = fe{2,8) = fc{3,4) = fc(4,5) = ^(5,9) = ^(6,7) = 1, all Other fc(,j) = 

(M, x') M = 6, x' = Je), e = 1, . . . , 6} = {(1, 9), (2, 8), (3, 4), (4, 5), (5, 9), (6, 7)} 



The crosslinkers are modeled as Gaussian springs, each connecting a pair of functional sites (i, j) (for i < j to avoid 
double counting). The network structure is uniquely specified by the list of N{N + l)/2 "link numbers" fc(ij) that 
specify the number of springs linking the pair of sites (i,j): x = {^(ij)} |14j . Then, the summation over network 
structures corresponds to summations over all (non-negative) integer values for fc(ij): 

oo CO oo 

E E ••• E - E ■ (2) 

X fe(i,i)=0 fc(i,2)=0 fe(N,]V)=0 
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Another, equivalent, specification of the network structure is given as foUows. First, we specify the total number 
of springs M . Second, we treat these M springs as if they are distinguishable, and number them using an integer 
index e = 1, . . . , Af. The e"^ spring can be denoted by an integer pair {ie,je) (with < j^) which specifies the pair 
of functional sites it links together. The network structure is then fully hut non-uniquely characterized by an integer 
M and a list of M integer pairs x' — {(*e, Je)}eii- Because all springs are identical, any arbitrary permutation of 
the M pairs within the list x' does not change the network structure. We easily find the number of distinct x' that 
correspond to the same network structure (which is uniquely determined by x), i-e., the "degeneracy factor": 

The summation over all network structures can then also be represented as a summation over the crosslinker number 
M, combined with a sum over the list of 2M integers {ie, je}, 

oo 

E=EE-- E • 

X M=Oll,jl IMJM 

Note, however, that the probability of a given network structure P^i must be divided by the degeneracy factor ([s]) in 
order to correct for the over-counting. We shall discuss this issue in detail below. For an illustration of a particular 
network structure, see Fig.[T] For the two equivalent specification of this network structure, see Table [Til 

Let us assert that after crosslinking the system is brought to the measurement state, in which the microscopic 
locations of the sites are now labeled by (rather than cv^ of the preparation ensemble) . The interaction between 
pairs of sites linked by springs is modeled via quadratic terms in the separations between sites: 

M 

^^xic) = 2p E l'^- - ''^f = 26^ E - '^^■l' • (4) 

This potential energy increases without bound as any two linked sites are separated, and hence it keeps them close 
to each other at all times. For later purposes, we also define the liquid-state partition function in the measurement 
state: 



Zuq = J i?ce-«'-(-). (5) 
The partition function and free energy of the crosslinked system are then given, respectively, by 



both of which depend on the network structure x- Throughout this paper, we shall set the Boltzmann constant 
to be unity. Note that, generically, _ffiiq(c) (i.e., the liquid Hamiltonian in the measurement state) differs from its 
preparation-state counterpart iJi'?q(c). We then average the free energy over realizations of the network structure x 
(with probability P^) to obtain 

^ = E,^x^x = -7^E,^xln^x- (7) 

Assuming that the selfaveraging property holds, this formula also gives the free energy of a typical sample. The 
calculation (or modeling) of the probability for the network structure x is the sine qua non for obtaining a proper 
statistical treatment of randomly crosslinked systems. This we now discuss in detail. 



III. GENERALIZED DEAM-EDWARDS DISTRIBUTION OF NETWORK STRUCTURE 



A. Statistics of connectivity between a pair of sites 



The crosslinking process that we shall consider is an instantaneous crosslinking process, as already mentioned 
previously. Many springs (crosslinkers) are introduced into the liquid system (in its preparation state), completely 
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FIG. 2: Two functional sites, Ci and Cj, each with an "efTective region" of radius e and volume Vt. Any end of a spring that is 
introduced into one such region is successfully linked to the site at the center of that region. 



randomly and independently of one another. We further assume that the extension of every spring d is distributed 
according to a Gaussian equilibrium Gibbs-Boltzmann factor with characteristic lengthscale b [cf. Eq. Q]: 

P(d)cxe-l'*l'/2^'. (8) 

In addition, we assume that each functional site has an "effective region" of radius e, centered around c". If an 
end of a spring is introduced at some position x that lies within the "effective region" around c^, it is successfully 
linked to this functional site. We assume that e is smaller than half the minimal distance between any two functional 
sites, so that any end of a spring cannot be linked to two sites simultaneously. This assumption does not limit the 
generality of our approach. 

Next, let pk (c°, c^) be the probability that a given pair of sites {c", c^} is linked by some integer number k of springs. 
If we assume that the crosslinkers are introduced into the system at random, independently and homogeneously, then 
k would be Poisson-distributed: 

p.(c^cO).C^e-h^sT/-^ (9) 

where we have introduced the dimensionless control parameter (m for the crosslinking density |15j , and the constant C 
is determined via overall normalization, see below. Note that Pkic^, c^) contains a factor proportional to the fc"^ power 
of Gibbs-Boltzmann factor ([s]), which characterizes the equilibrium distribution of each spring. It is convenient to 
define a (three-dimensional) soft delta function db{x) in terms of a Gaussian spatial profile having variance parameter 

Sb{x) = {V2^/bfe-^"/^''\ (10) 

The function Sb{x) tend to the Dirac delta function as 6 — > 0. The linking probability ([9| can then also be expressed 
in the following concise form: 

where we have exchanged the crosslink density control parameter fi for the parameter /i, defined via 

^={b/V2^)^ fi. (12) 



It is easy to check that Eq. (11) is properly normalized: J2^^o Pk{Ci , c'j) = 1. 



B. Statistics of the connectivity of the network 

We assume that the liquid, prior to crosslinking, is in thermal equilibrium with respect to Hamiltonian H^^^{c'^), 

where c" = {cf} denotes the locations of all N functional sites in the preparation state [16]. The joint probability 
density function for all N sites (prior to crosslinking) to be at the locations {cf, . . . , c^} is then given by 

P(c°) = ^e-<(^"), (13) 

^liq 
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where the normaUzing partition function is given by 

Z^,^ = I Dc''e'""^^^"\ (14) 

Now recall that the network structure is characterized by a set of A^(A^+l)/2 stochastic variables [T7|, or link numbers, 
j)}, each of them characterizing the number of springs linking a particular pair (i, j). Given a liquid configuration 
prior to crosslinking, these link numbers are fully uncorrelated, and governed by the following conditional probability 
distribution: 



n (A^,(cO - ^))''-'^ , (15) 



where we have used the linking probability given by Eq. ( |11[ ). Making the further definition 

-ffnorm(c) = -/2 ^ 5fc(Ci - Cj) , (16) 

we can write this conditional probability as 

P(x|c") = (-")nr^(A^.(c?-c°))'=<-' 



^— ' n kT-f 

ff„o™(c")-Aff,(c'') TT A'<-' 



e 



n 



(1,1) ■'■•'> 



where we have used Eq. (10 1, Eq. (12 1, and Eq. ([4| in the above three equalities, respectively. If we were to specify 
the network structure using the set = {(*ei Je)j e = 1, • ■ ■ i M} instead of the set x, the corresponding conditional 
probability would be 



P(X'|C°) = (-")^ n ^'5^(4 - 4)' (18) 



e=l 



where we have divided the probability by the proper degeneracy factor ([3|, so as to cancel the over-counting. 
The probability of obtaining the network structure x can be arrived at by using the law of total probability as 

P^ = ^P(x|cO)F(c°) 
= i / i?c"e-^".(^")+«— (^'')-^^>^(^'') TT ^ . (19) 
Furthermore, by making the definition 



we obtain 



Hence we see that, apart from a dimensionless factor, is the ratio of two partition functions. Even though the 
variables (associated with preparation state) appear only as dummy variables, the integrations over them in Eq. ( 20 ) 
cannot be carried out explicitly. Below, when we calculate the disorder average of free energy, we shall find that these 
variables interact with the corresponding variables in the measurement ensemble. It is this interaction that 
generates the ubiquitous memory effects in randomly crosslinked materials. 
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C. Normalization of probabilities 



The conditional probability that there are precisely M linking springs in the system can be obtained by summing 
Eq. (18) over all possible values taken by the indices {ie,je)' 

N N N 

2l,il = l 12 j'l=2 iM.jM=l 

= ^"--'''"'^^-H^oUc'))'' ■ (22) 
It is elementary to check that the conditional probability P{M\cq) is properly normalized for any fixed Cq: 

oo 

E P{M\c^) = 1. (23) 

M=0 

To check the normalization of the probability P^, we sum Eq. ( [l9|) over all possible network structures. By Eq. ([2|, 
this amounts to summing over all possible linking numbers {^(ij) }• Further using Eqs. Bl we arrive at 



^Uq^ X (7,7) ^('-J')' 



The summation and product can readily be calculated, leading to the factor exp — i7x(c*'), which precisely cancels 
the corresponding exponential factor. Hence, as expected, we have the normalization condition 

E^x-1- (25) 



D. Averaging over network structures 

We are now ready to compute the disorder-averaged free energy, Eq. Q. For this purpose, we use the following 
identity, which is commonly called the replica trick |18j : 



\nX = lim - InX". 

n-i-O n 



(26) 



We use this to rewrite the average free energy as 

F = -T lim - InZ^ = -T lim - In V P^Z" 



-T lim - In y n 



(27) 



We calculate the right hand side of Eq. (271 for positive integers n, and subsequently analytically continue n to real 



values. More specifically, we calculate ZI^ by replicating the variables {ci} a number n times, to obtain an n-tuple 
,c"). Together with the corresponding variable of the preparation ensemble, they form a (1 + n)-tuple 
,c"). The shorthand c has been extensively used for this (1 + n)-tuple in the vulcanization theory 



(c°, 



literature. Using Eqs. (|6 20 1, the product of partition functions Zq^^{Z^)"' in Eq. (271 can be written in terms of 



functional integrals (with tHe shorthand notation Dc = Hq^o -De") as follows: 



^0,x (^x 



DC XI 



(28) 



c«=0 
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Next, by using the definition Eq. Q we easily see that 



a—Q a— (.^J) 

^ Ai/x(c), 



(29) 



where we have used the notation |c"|^ = X]a=o I'^f 1^- "^^^ ^'^^ substitute Eq. (28) back into Eq. (27), exchange 
the order in which and J Dc are performed in Eq. (27), and sum over x foi' ^ fixed configuration of c. The 
summation over the reahzations of network structure can be calculated using the following result: 



n 



E n 



= nE 

[id) k 
(id) 

= exp^^ ^ Sb{ci - Cj) 
where in the last equality, we have used the definitions 

71 

5{x) = l[S,{x'')^5,{x")Sb{x^)---6b{x"). 



(30) 

(31) 

(32) 
(33) 



Q = 



As we shall eventually take the replica limit, n — ^ 0, the difference between fj. and p, [see Eq. (12)] can be ignored. 



Combining all these results, we find that the disorder-averaged free energy is given as follows: 

F 



-T lim - In ( -jj- Zi+„ ) , 



7?i+„(c) - i/Sq(c°) - i/„orm(c'') + ^ i?liq(c") + i/„oon(c) 



N 



(id) 



(34a) 
(34b) 
(34c) 

(34d) 



The Hamiltonian -ffnorm[c] represents an effective short-range attraction between all pairs of functional sites in the 
replicated coordinate space c = {c"}^^oi resulting from cross-linking. This short-range interaction reduces to a delta 
bmction in the limit 6 — >■ 0. 
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IV. PHYSICAL QUANTITIES AND CORRELATION FUNCTIONS 

We shall use the notation (^(c))i+„ to indicate averages of physical quantities with respect to the weight in the 



replicated partition function Eq. ( 34b I 



= J Dee- " J • 

When there is no risk of confusion, we shall neglect the subscript 1 + n on the average and simply write (^'(c)). It is 
understood that at the end of the calculation of such quantities the replica limit n — should always be taken. 



In the replicated Hamiltonian Eq. (34c I, only the last term, 77noim(c), couples together various replicas. It 

can be understood as a short-range attraction between all pairs of the 1 + n replicated positions c = (c°, c^, . . . , c"). 
The strength of this attraction is controlled by the density of crosslinking fj,. If = 0, i.e., the system is not crosslinked, 
HnoYui — Oj and reduces to a sum of liquid Hamiltonians, one for each replica: 

n 

Zl+n ^Uq ■ (-^liq)" ; 

F ^ -TlnZiiq. 

Variables in distinct replicas then become completely decoupled from one another, and the replicated partition function 
Zi^n describes the equilibrium liquid in the preparation state together with n copies of equilibrium liquid in the 
measurement state, with no interactions amongst these liquids. 

For systems with crosslinks, /i 7^ and so there are nonvanishing correlations between the variables c° and the other 
c". These correlations characterize the dependence of the properties of these systems on the history of formation. Let 
us look systematically at the various physical quantities that can be constructed within the framework of an (1 + n)- 
replica theory, and relate them to the suite of correlators involving the preparation and measurement ensembles 
discussed in the Introduction. Below, we use the shorthand notation A" = A{c") for physical quantities, and we 
frequently omit the subscript 1 + 71 on various expectation values, all to ease the notation. 

1. One-replica quantities: {A*^) and (A") {a 7^ 0) are average physical quantities in the preparation state and 
in the measurement state, respectively. In the notation used in Sec. jljthey are written as [A]p and (A)r 
respectively. 



There are also correlation functions in a single replica: (A'^ = [A-B]p and {A" B") — (AB) 
which quantify correlations in the preparation state and in the measurement state, respectively. 



[a ^ 0), 

p 



2. Two-replica quantities: {A" B^) — {A)n-i{B)„-^ (for a ^ f3;a, /3 0) are nonvanishing in the presence of the 
quenched disorder arising from random crosslinking, and thus give the glassy correlations in the measurement 
state. The order parameter for the random solid state is an especially important example of such a multi-replica 
quantity. 



Memory correlation functions involving 0-th replica only once, {A^B")c = A (-B)m (for a ^ 0) measure the 

P 

cross correlation between the preparation state and the measurement state. These memory effects are a hallmark 
property of randomly crosslinked materials. Vulcanization theory provide the unique framework via which they 
can be systematically calculated. 

In the setting of liquid crystalline elastomers, the memory correlators are the quantities that can best distinguish 
systems prepared under distinct physical conditions, e.g., isotropic genesis nematic elastomers (IGNEs) and nematic 
genesis nematic elastomers (NGNEs). Various correlation functions, expressed in both replica notation and physical 
notation, are summarized in Table IIVI 



A. Significance of the zeroth replica 

We can use Zi+n with appropriate source terms to calculate the average of a physical quantity A(c°) (= A°) that 
only depends on the degrees of freedom c" of the preparation ensemble. As replicas 1, . . . ,n play no role in such a 
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TABLE II: Various correlation functions in vulcanization theory, showing the correspondence between physical notation and 
replica notation. (Note that in the table we have 7^ a 7^ /3 7^ 0.) In the replica notation, all averages are defined by Eq. (35 1. 



Preparation Measurement Glassy Memory 

correlation correlation correlation correlation 



Physical 
notation 



notation 



{A^B") {A"B°') 



calculation, we can take the replica limit n — > before the calculation of the expectation value. This is equivalent to 



deleting all variables c" for a = 1, . . . ,n from the replicated Hamiltonian (34c I. In particular, we have 



-^noi'm(^) ^ -^norm('^ ); 

i7i+„(c) ^ i?" (c°). 



We also have to delete the integrals over c" (for a = 1, . . . ,n) in the functional integrals (34b) and (35). Therefore, 
we have 







As a result, {A°) is indeed the average in the preparation state, as expected. 



(36) 



B. Replica limit and causality 



In the replicated partition function, Eq. (34b), there is one copy of the preparation ensemble (the 0-th replica), 
and n copies of the measurement ensemble. In Fig. [3] we show a cartoon that shows schematically the relationship 
between various replicas. All I + n replicas interact with one another through the last term in Eq. (34c), which is 
linear in the crosslinking density. Note that there are n copies of the measurement ensemble (i.e., replicas 1, . . . ,n) 
but only one copy of the preparation ensemble (i.e., replica 0). In general, owing to the interactions, the preparation 
and measurement ensembles can mutually influence e ach other. For example, if we tune some parameter J° in fluid 
Hamiltonian iJf?^ of the preparation state [see Eq. (34c|], there will be changes in the physical quantities in the 
measurement state: 



d 



(^)m ^ 0. 



(37) 



This is, of course, entirely reasonable, as the properties of polymer networks do depend on their history of formation. 
Reciprocally, however, if we change a parameter J in the measurement state, might there, at least in principle, also 
be responses of physical quantities in the preparation state, e.g. 



d_ 
dJ 



■A],^0 (?). 



(38) 



Such a response would be entirely unphysical, as it would violate the principle of causality. In the framework of 
vulcanization theory, the causality principle is rescued via the replica trick: as n — 0, the weight of the measurement 
ensemble relative to the preparation ensemble tends to zero, and therefore the former can have no quantitative 
influence on the latter. 

An explicit proof of this result is simple. Let the parameter J couple linearly to some quantity ^(c) in the (1 + n)- 
replica Hamiltonian i?i+„; this leads to the change 



a=l 
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FIG. 3: Cartoon showing the logical relations between different replicas: the preparation ensemble and n copies of the mea- 
surement ensemble are coupled to one another via the crosslinking Hamiltonian //norm(c). The relative weight of measurement 
ensemble vanishes in the replica limit, thus protecting the causality principle. 



Substituting this back into the expectation value Eq. 
J — 0, we obtain 



( 36 ) , taking the derivative with respect to J, and finally setting 



(39) 



where in the last equality we have assumed that permutation symmetry amongst replicas 1 to n remains intact. The 
right hand side of Eq. (39 1 therefore vanishes in the replica limit, i.e., n ^ 0. Even in the presence of the spontaneous 



breaking of this replica permutation symmetry, the right hand side of Eq. ( 39 ) would even be proportional to n, and 
hence would vanish in the replica limit. 



Deam-Edwards revisited 



In the replicated Hamiltonian (34c), the parameters of the preparation state only appear in H^^^, whereas the 
parameters of the measurement state only appear in iJ^q. These two types of parameters can be separately controlled 
in experiments, but they jointly determine the statistical physics of randomly crosslinked materials. Let us consider 
a rather special case in which i/f?^ = Hnq, i.e., a system that is prepared and measured under identical conditions. 

What appears inside the (1 + n)-replica Hamiltonian (34c) is, however, iJj"q(c°) — i/norm(c") instead of just iJj"q(c"): 



Hl+„(C) = Fnq(c°) - H„orm(c°) + J2 ^liq(c") + H„orm(c). (40) 

a=l 



Consequently, the replica Hamiltonian Hi^„ is not invariant under permutations of the replicas that mix the 0-th and 
the other replicas. Hence, the average of a quantity A in the preparation state is generically different from its average 
in the measurement state: 



This is completely reasonable, as the random network structure does affect physical quantities in the measurement 
state but cannot not influence the preparation state. One simple example involves the nematic correlations in liquid 
crystalline elastomers that have been crosslinked in isotropic state. The random polymer network tends to disorder 
the nematic degrees of freedom, and therefore reduces the nematic correlations. 

Historically, Deam and Edwards [T] chose a particular probability distribution for the network structure given 

by 



pDE ^ fo^ TT \^ UX) 

liq (j J) 



where 



(42) 
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The corresponding (1 + n)-replica Hamiltonian is then given by 

n 

Hff„(c) = i7liq(c") + 5]i?liq(c")+i7„orm(c), 
n 

= ^i/liq(c")+ffnorm(c). (43) 

Q = 

in which all 1 + n replicas appear symmetrically. The resulting 1 + n replica theory then has an enlarged 
permutation symmetry, rather than the Sn permutation symmetry that is mandatory for replica theories. The Deam- 
Edwards choice of disorder distribution is appealing, in that it simplifies the analysis considerably. However, it also 
leads to the unnatural consequence that the averages of arbitrary physical quantities in the preparation state are equal 
to their counterparts in the measurement state. The S'l+n permutation in the Deam-Edwards theory is therefore not 
mandatory, as first pointed out by Broderix et al. [4^. The explicit breaking of the 5'i+„ permutation symmetry down 
to Sn permutation symmetry is physically natural, and opens the way to a systematic analysis of the interrelationship 
between randomly crosslinked materials and their histories of formation. This is a valuable and important element of 
the physics contained in vulcanization theory. 



D. Comparison with spin glasses and other disordered systems 



We end this section with a brief discussion of the difference between the vulcanization model and certain well-known 
spin glass models. Consider, e.g., the Edwards- Anderson model W: 

H = — ^ JijSiSj. 

The quenched random variables {Jij} are typically assumed to have independent Gaussian statistics with vanishing 
means, and therefore averaging over them (after application of the replica method) is elementary. The final, replicated 
theory reads: 

(id) a,p 

and contains n replicas of the annealed degrees of freedom (i.e., the spins). The variance of the quenched variables tJj 
appears only as a parameter in the replicated Hamiltonian. The quenched variables {Jij} have been integrated out, 
and therefore do not show up in the replica theory. By contrast, in vulcanization it is not just expedient but indeed 
manifestly physical to employ an equilibrium preparation state to generate the statistics for the network structure. 
(For vulcanization, the quenched random variables are the network structure.) Furthermore, correlations between the 
preparation state variables and the measurement state variables encode the permanent memory effects, and these are 
experimentally measurable. Finally, we note that the Edwards- Anderson model represents a considerable idealization 
of disordered systems in real world. In general, the physics of most disordered systems does depend on the history of 
formation. Hence, the proper statistical modeling of such systems requires a separate statistical ensemble. 



V. ORDER PARAMETER FOR VULCANIZATION 



The crosslinking-induced interaction (34d) couples distinct functional sites. It can be decoupled using the standard 



Hubbard-Stratonovich transformation. To proceed, we first notice that Gaussians are closed under convolution. 
Therefore, the following identity about the (1 -|-n)(i-dimensional soft delta function Sb{x) [defined in Eqs. ( 10|33 )] can 
readily be established: 

5b{x-y)= J dzSb'{x - z) Sb'iz - y), b'^b/\/2. 
Equation (34d) can then be written as 

» TV N 
Hnoimic) = dz^ 5b' [Ci ^ 5b' {Cj - z) 

— — ^ J diw(£,c)^, 
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where 

N 

uj{z,c) = ^6t' {c,- z) (44) 

is the microscopic collective density of functional sites, coarse grained by the Gaussian wave packet db'{x). It depends 
explicitly on the microscopic configuration of the replicated system. Our definition of the collective field differs from 
those in the other literature on vulcanization theory [5] via a multiplicative factor of N. We can now introduce an 
auxiliary field il,{x) = il{x'^, . . . , a;") via the following Hubbard-Stratonovich transformation: 



where C is a trivial constant that will be neglected henceforth. Insertion of this result into Eq. (34b) leads to 

Z^^^_^ = j Dc j i:)f7(£)e--f^nd['='Vffx[c«]-ES=l^nJc°]-M/dxO(*)^+2A./dxa;(x)f2(i) 

= j i:»f7(A)e"^-'''^*'^^(*)'y" i:)ce~-f^Hi['=°i+-f^x['=''i-s:s=i^^u,[c°]+2M/'i*"(*)f2(*), (45) 

Let us further define a modified liquid state partition function via 

= j Dc'' £-^".['=''1+^- , (46) 

and the average ( • • • )5'+„ as 

(*(c))?+„ ^ ]^ / i^c*(c)e-«-[^"l+«-[^"l-E:..«,TJc"l. (47) 
^liq ■ (.'^liq) J 



Then we can rewrite Eq. ( 45 ) as 

= Zl^ ■ (Ziiq)" • j Dn{x) e-^-("), (48) 
where the effective Hamiltonian iJoff (f^) is defined as 

i7^ff(f]) =^y"d£17(A)2-ln^e2^-/''^*'^(*)^^(*)^"^ . (49) 
We can now define the average of functions of il using the equilibrium partition function associated with Hcsi^)'- 

The following identity is a well-known property of the Hubbard-Stratonovich transformation: 

(f7(A)),f, = (^i,c))i+„ = /f^5fc,(^°-c°)---<5,,(a;"-cr)\ , (51) 

\i=l / 1+n 

which shows the physical significance of the collective field r2(A), namely that its average is the joint probability density 
function of the following 1 -I- n events: in the preparation state there is a particle at some point a;", whereas in n 
independent measurements on the measurement state the same particle is found at positions a;^, . . . ,a;", respectively. 
Using the physical notation discussed in the Introduction, this probability density can be expressed as 

N 

("(*)>cff = E - c,)(5fc'(a3i - c,)),„ • • • (<5b'(a;" - c,))„J , (52) 

where, as usual, the subscripts p and m respectively denote the preparation and measurement states, and the overline 
indicates an average over realizations of crosslinkers. U,(x) encodes information about the localization of particles in 
the gel phase, and consequently is called the vulcanization order parameter. The definition of ^l{x) in this work differs 
from that in some of the other vulcanization literature via a trivial additive constant and a multiplicative constant 
N. 
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A. One replica sector and marginal distributions 



We can integrate out the n variables x^,. 
thereby obtam a one-rephea distribution 



x" in the order parameter (51), leaving only one variable a;'', and 



{n"ix^)),ft 



„ n IN 

/ \{dx'^{n{x))^lS25A^'~^.) 

a=l \j=l 



(53) 



l+ri 



The resulting object is the average number-density of functional sites in the preparation state, which is manifestly an 
intensive quantity. At the coarse-grained level, it can also be understood as the average particle density of the liquid 
prior to crosslinking. The fluctuation 5Vf'{x'^) — ^^{x'^) — (il"(a:"))cff is also linearly related to the fluctuation of 
particle density in the preparation state. Similarly, we can obtain the a-th replica distribution for a 7^ 0: 



{n^{x^)) 



off 



[| dxf^{Vl{x)) 




(54) 



which describes average particle density in the measurement state. This relationship between the vulcanization order 
parameter and the particle densities in the preparation and measurement states is the primary reason that we choose 
the particular normalization of the order parameter field in Eq. 



(44). 



For swollen gels and elastomers there can be large density fluctuations in the system (either quenched or thermal). 
These fluctuations, characterized by correlations of the one-replica distributions il"(a;"), may strongly influence 
the elasticity of system. For unswollen gels and elastomers, the mean particle density at long enough lengthscales 
is usually uniform and its fluctuations are negligibly small. The one-replica part of free energy then contains no 
interesting physics, besides the stabilizing of the density fluctuations. 

Finally, by integrating out n — 1 variables, we obtain the two-replica correlator 



{n''^{x"-xP)),s^ 



n dx^nix)), 



(55) 



which encodes the system average of the correlation of the position fluctuations of a site in two independent measure- 
ments. Macroscopic translational symmetry dictates that this function depend only on the difference between two 
coordinates, x" and x^ . Generically, it contains less information than the full order parameter Q,{x). 



VI. LANDAU THEORY AND SADDLE-POINT APPROXIMATION 



Eqs. (48 [50 ) provide a formal approach to calculate all relevant physical quantities. Our problem is therefore 
reduced to the calculation of a functional integral over 17 (a;). At the saddle-point level we approximate the disorder- 
averaged free energy ( 34a ) using the minimum of the effective Hamiltonian i?off [^] , Eq. ( 49 ) . Even after making this 



approximation, the problem remains difficult, because we do not know the effective Hamiltonian (49) exactly, let alone 



its minimum. We therefore proceed to e xpa nd the exponential in Eq. (49) in a Taylor series in il. At this stage, it is 
convenient to shift the order parameter (51 ) by a constant N/V^^" (where TV is the total number of functional sites), 
so that the it satisfies the condition 



J dxfl{x) 



0. 



(56) 



Both in the liquid phase and in the gel phase, the system is uniform in density. Therefore the saddle-point order 
parameter satisfie: 



Yl dx^n{x) = [(p(a=")). 



N 
V 



0. 



(57) 



We therefore only have to consider the subspace of fl{x) in which Eq. (57) is obeyed. This subspace is usually called 
the higher replica sector. 



The form of Landau free energy given in Eq. ( 49 ) is generic for all isotropic networks 



Hhr — 



I (1=0 



-rn^ 

2 



3 



(58) 
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The coefRcients in the Taylor expansion do of course depend on the short-scale structure of the constituent particles. 



It may appear surprising that the cubic term in Eq. (58) has a negative coefficient, which seems to suggest that the 
Landau theory is unstable, at lest for uniform values of the order parameter. By Eqs. ( 53|54| ), however, a large and 
uniform n{x) also implies the large deviation of the particle densities from their equilibrium values in each replica [and 
therefore would violate Eq. ( 57 1] , and this would incur a large free-energy penalty. The Landau theory of vulcanization 
is therefore locally stable. Minimization of Eq. ( 58 ) leads to the following saddle-point equation: 



E 



(59) 



This equation should be solved subject to the constraint (57), and is therefore necessarily nonuniform in replicated 
space. 



A. Order parameter at the saddle-point 



By solving the saddle-point equation (59), subject to the constraint ( [57[ ), we obtain the saddle-point value of the 
order parameter, which we denote by ^l(x). For r > 0, the saddle-point value of fl is trivially zero, corresponding to 
the liquid state. For r < 0, the trivial saddle-point becomes unstable, and a nontrivial (i.e., replica-space dependent) 
saddle-point value emerges, corresponding to the amorphous solid state. It has the following form of superposed 
Gaussians [T^: 

n{x) = g JdzJ dTp(T) (2vrr)(^+")''/'e-5 ^==o(^°-^)' ^. (60) 

The coefficient q is identified with the number density of the infinite cluster (i.e., the gel fraction). The dummy 
parameter r is the inverse variance of the Gaussian fluctuations from their equilibrium positions of the localized 
particles (belonging to the infinite cluster), and is called the inverse square localization length. The fact that there is a 
distribution of r signifies the heterogeneous nature of randomly crosslinked systems; the infinite cluster is more rigid 
in some places and looser in others. Both q and p{t) are determined by solving the saddle-point equation. For details, 
see Refs. [BI[7|. It is important to note that they are independent of strain deformation and all other parameters that 
are tunable in the measurement state. The physics of the saddle-point order parameter can be better understood by 



focusing on the integrand inside Eq. ( 60 ) : The vector z can be interpreted as the mean position of a particle in the gel 
fraction, and is localized around x'^, the position of the particle at the moment of crosslinking. After crosslinking, the 
particle in the gel fraction (x", a = 1, . . . , n) is localized around z, its mean location in the measurement state. The 
uniform integral over z implies that the gel is statistically homogeneous, after averaging over the quenched disorder. 

If the crosslinker density is not high enough, the system is in the sol phase in the measurement state, no infinite 
cluster emerges, and all particles are delocalized, so that there is no correlation between their positions in different 
measurements, performed in the preparation and measurement states. Hence, the associated joint pdf is a trivial 
constant, and the average order parameter is vanishes identically. By contrast, in the gel phase, there is an infinite 
cluster which constitutes a finite percentage of the overall mass. For any particle in this infinite cluster, different 
measurements of its position are necessarily correlated: knowledge of its location in the preparation state tells us 



information of its whereabout in the measurement state, and hence the average order parameter in Eq. (51) is 
nonzero. The order parameter ft{x) thus distinguishes the gel phase from the sol phase. 

We can actually deduce the functional form of the saddle-point order parameter ([60| using the definition of the 



order parameter (52) [after subtracting a trivial constant, so that the constraint (56) is satisfied] together with some 
simple, intuitive reasoning. We have already noted that particles belonging to the iquid fraction do not contribute 
to the saddle-point order parameter. For a particle i belonging to the infinite cluster, we expected that its position 
in the measurement state is Gaussian-distributed around its mean value, which we denote by Zi, with variance l/r^. 



The thermal averages inside Eq. ( 52 ) combine to give 



{St-ix' - c,)),„ • • • {Sb'ix^ - c,))^ = (2^T,)"'^/2 exp 1^" - ^'I'l • (61) 

The parameters Zi and of course depend on the microscopic configuration of the preparation state at the instant 
of crosslinking, as well as on the crosslinks that are realized. It is important to note that the mean position Zi is 
distinct from the position of the particle at the instant of crosslinking, which we denote by c^. Although, for a given 
configuration of preparation state and crosslinker realization, Zi is fully determined, it becomes a statistical variable 
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after we include the various crosslinker realizations. In fact, we expect that Zi is Gaussian-distributed around the 
particle position in the preparation state, with the same variance l/r^. Also, for different realizations of crosslinkers, 
the variance l/r^ itself should also be different, obeying some distribution p(t). An ergodic hypothesis suggests that 
this distribution is reflected in the microscopic spatial heterogeneity that individual samples of elastomer exhibit. All 
these considerations suggest that 



{bv{x^~c,))^---{by(x^~cC))^ = Jdz J dr(27rr)(i+")'^/2p(T)exp|-^f^^ 



(62) 



Finally, we note that in the preparation state the system is a homogeneous liquid and thus translation invariant. 
Hence, the position of a given particle Cq is uniformly distributed over space, and therefore for an arbitrary function 
/(c°) we have 



(63) 



We now substitute Eq. ( 62 ) back into Eq. ( 52 ) and carry out the average over the preparation state and also sum over 



all of the particles. Thus, we finally arrive at the precise form of the saddle-point order parameter given in Eq. (60 1 



B. Translational symmetry and order parameter 



The 1 + n replicated Hamiltonian Eq. (34c I, possesses the symmetry of independent translations of the 



replicas, i.e., it is invariant under the a-dependent translations: 

cf^< + u", i^l,...,N. 

This symmetry is also possessed by the average order parameter n,{x) in the liquid phase (which has the value zero). In 
the gel phase, however, this translational symmetry is explicitly broken by the saddle-point order parameter, Eq. ( |60[ ). 
Translations of all particles in one replica leads to a distinct but energetically degenerate order parameter: 

n{x°, ...^x" +u,... x") ^ n{x°, . . . , a;", . . . x"). 

This reduction in in symmetry is a consequence of the localization of the particles associated with the emergence of 



an infinite cluster. Nevertheless, translations of all replicas by a common vector u remains a symmetry for Eq. (60 1: 

n{x°, x\...,x"-)= n{x° +U,X^ +U,...,C'' +U). 

This macroscopic translational invariance (MTI) ^ reflects the fact that gels and elastomers are statistically ho- 
mogeneous. This symmetry can be broken by translational ordering, such as smectic ordering in liquid crystalline 
elastomers. 



C. Classical theory of the elasticity of rubber 



Because rubbery materials can typically sustain large deformations, their elasticity theory is necessarily nonlinear. 
This is already apparent in the classical theory of rubber elasticity, which gives the elastic free energy 

i7(A) = ^TrA^A, (64) 

in terms of the (spatially uniform) deformation gradient matrix A. The matrix A relates the undeformed positions 
of particles belonging to the infinite cluster to their deformed positions z' = A ■ z. As most rubbery materials are 
nearly incompressible, det A can be taken to be unity. This is largely responsible for the nonlinear nature of rubber 
elasticity. In the framework of statistical physics, the positions of particles fluctuate, and this deformation relation 
must be understood in the sense of relocation of the mean positions. It is reasonable to suppose that the saddle-point 
order parameter of the deformed state is given by 

n^{x) ^qJJ^ 1-°-^^!^ - (65) 
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Note that the mean position z is afBnely deformed in the measurement state (repHcas l,...,n), but not in the 
preparation state (rephca 0), and that q and pij) are not deformed (i.e., do not depend on A). Deformation of 0-th 
rephca can always be "gauged away" via a coordinate transform of the dummy integration variable 2, and hence has 
no physical significance. It is the relative deformation between the preparation and measurement ensembles that has 
physical significance. 



It has been shown explicitly [7] that Eq. ( 65 ) indeed satisfies the saddle-point equation and the boundary conditions 
for a solid having a uniform deformation A. Furthermore, by inserting Eq. (65) into the Landau free energy ( |58[ ) [ 20] 
it is straightforward to obtain the elastic free energy of the deformed system, as given in Eq. (64), i.e., the classical 
theory of rubber elasticity [5|. Moreover, one finds that yit cx |r|^ for the shear modulus [Hlj. It is rather interesting to 
see that classical rubber elasticity theory emerges at the saddle-point level of vulcanization theory. 



D. Nematic elastomers and neo-classical elasticity theory 

In nematic elastomers, the polymer chains carry liquid crystalline mesogens that are prone to undergoing nematic 
ordering. The mesogens can either be parts of the polymer backbone (i.e., main-chain nematic polymers) or attached 
to the backbone side-on (i.e., side-chain nematic polymer). A Landau theory of nematic elastomers has be derived 
via vulcanization theory elsewhere [7]; here we simply quote the results. 

Let the (spatially uniform) nematic order parameter be in the preparation state and Q in the measurement 
state. The free energy ( [58| needs to be modified to include possible couplings between the nematic and vulcanization 
order parameters. To lowest order in these order parameters, the couplings are 



(66) 



These couplings can be incorporated into the isotropic term in Eq. ( 58 ) to yield 



Hi? 



1 I yof] + ^ |] vjn + ir - + . . . | , 



where 



f = I + Q" and I 



Q 



(67) 



(68) 



are the polymer step-length tensors in the preparation and measurement states, respectively. 

The general saddle-point in the presence of a macroscopic deformation gradient A, which we denote by Qa.{x), can 
also be obtained, by minimizing the free energy (67): 



(69) 



where 



= x° -z and yl = x°' -A- z. 



(70) 



Note that the Gaussian fluctuations of the particle positions have variance matrices in the preparation state and I 
in the measurement state. The effect of nematic order is therefore to render the thermal position fluctuations of each 
molecule anisotropic. 



Substituting this deformed saddle-point (69) back into Eq. (67), and carrying out the prescription specified in 

(71) 



Eq. (34a), we find the following elastic free energy for nematic elastomers: 

^ Tr Z"ATr lA, 



which is precisely the neo-classical elasticity theory for nematic elastomers, derived by Warner and Terentjev [31 15] as 
a generalization of the classical theory of rubber elasticity. 
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VII. CONCLUDING REMARKS: WHY VULCANIZATION THEORY? 

It is quite satisfactory that the classical theory of rubber elasticity and the neo-classical elasticity theory of ne- 
matic elastomers can be derived from the vulcanization theory as saddle-point approximations. These theories were 
originally established by studying the statistics of single polymers inside a network, with two key assumptions [31 [5] : 
(1) that strain deformations are affine at all scales; and (2) that the polymer statistics are Gaussian. These theories 
have attained classical status because of their remarkable successes in explaining the basic properties of isotropic 
and nematic elastomers. Nevertheless, the starting point of single polymer statistics substantially restricts the gen- 
eralizability of these theories. For example, it is hard to see how one might systematically explore the nature of 
heterogeneities in random polymer networks using these molecular-level approaches. 

Vulcanization theory was developed to understand the heterogeneities of randomly crosslinked systems. In the face 
of its (perhaps heavy!) machinery — such as field theory, the replica technique, and multiple statistical ensembles — if all 
that could be obtained from it were re-derivations of classical results then it would be hard to argue that vulcanization 
theory has been worth investing in. Fortunately, however, several fundamental new results have already been obtained 
using vulcanization theory, concerning the heterogeneity of rubbery materials: (1) the distribution of localization 
lengths p{t) in isotropic networks [lU]; (2) the distribution of internal random stresses in isotropic networks pH, I12| : 
and (3) the statistics of random fields and memory effects in isotropic-genesis nematic elastomers [T3]. We refer the 
reader to the literature for further details about these topics. We think it is fair to say that we are still in an early 
stage in the unraveling of the statistical physics of network heterogeneities and memory effects in randomly crosslinked 
materials. We hope that fresh and interesting results will be available to report in the near future. 
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